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Abstract 

We consider a two dimensional lattice model to describe the opening of a crack 
in hydraulic fracturing. In particular we consider that the material only breaks 
under tension and the fluid has no pressure drop inside the crack. For the case in 
which the material is completely homogeneous (no disorder) we present results for 
pressure and elastic energy as a function of time and compare our findings with 
some analytic results from continuum fracture mechanics. Then we investigate 
fracture processes in strongly heterogeneous cohesive environments. We determine 
the cummulative probability distribution for breaking events of a given energetical 
magnitude (acoustic emission). Further we estimate the probabilty distribution of 
emission free time intervals. Finally we determine the fractal dimension(s) of the 
cracks. 
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1. Introduction 

Hydraulic fracturing is widely used in soil mechanics to improve the perme- 
ability of reservoirs either in oil recovery or of geothermal wells An incompress- 
ible fluid, in general water, is pushed under high pressure deep into the ground 
by injecting it into a borehole. The fluid penetrates into the solid opening long 
cracks. In order to optimize this rather costly procedure it is crucial to get some 
deeper understanding of how the fracturing occurs. 

In the fleld it is unfortunately very difficult to obtain direct information of 
the crack evolution in the ground. In present engineering essentially two types 
of measurements can be performed: On one hand one can monitor the pressure 
fluctuations at the injection pump and on the other hand one can register acoustic 
emission signals on the surface. In two dimensional Hele-Shaw cells some controlled 
laboratory experiments have been performed by injecting water or air into the 
center of the cell. The resulting cracks display a ramifled structure which for high 
enough pressures is fractal with a dimension of 1.4 - 1.5. 

Using a triangular network of springs and radially stretching the network 
on the outer boundary into the six directions of a hexagon the breaking of a 
material from a central hole was investigated by several authors . They observed 
fractal cracks having a dimension that depended very much on the type of applied 
displacements (shear, uni-axial, radial). A stability analysis t'^^ of the boundary 
of a circular hole with internal pressure has, however, shown that this case differs 
considerably from that of a stretched membrane due to the non-linear dependence 
of the growth velocity of the crack surface arising from the threshold in cohesion 
force that must be overcome to break the material. 

We have introduced a model in which the imposed load represents a pres- 
sure that acts along the entire (inner) surface of the crack in a direction perpen- 
dicular to the surface (von Neumann boundary value problem). In this way, the 
point of application of the imposed load varies during the growth of the crack, 
a situation that describes the case of hydraulic fracturing more realistically than 
previous spring models. It also turned out to be more efficient to use a beam 
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model instead of springs. This model had, however, two drawbacks: On one hand 
the pressure was kept fixed while in real applications it is usually easier to sus- 
tain a fixed injection rate. On the other hand the heterogeneities of the medium 
were "annealed" , i.e. changing in time while the disorder in breaking strength or 
stiffness in real soils is usually "quenched" , i.e. constant on the time scales of the 
breaking process. 

In the present paper we present a model with constant injection rate in which 
the variations of pressure can be measured and in which the cohesion force is a time 
independent random variable. We investigate the strong pressure fiuctuations and 
measure the energy release as function of the statistical distribution of cohesion 
forces. 

2. The Model 

In the following we will outline the employed model. First we give a brief 
description of the basic elastic equations and an explanation of how to incorporate 
heterogeneous cohesion properties into the fracture model. After this we explain 
in detail the used boundary conditions. Finally we present the employed breaking 
rules which contain the physics of the considered breaking process. 

We consider the beam model (as defined in p. 232 of Ref. 6) on a two dimen- 
sional square lattice of linear size L. Each of the lattice sites i carries three real 
variables: the two translational displacements Xi and j/, and a rotational angle 
(fi. Neighbouring sites are rigidly connected by elastic beams of length I. The 
beams all have the same cross section and the same elastic behavior, governed 
by three material dependent constants o — 1/{EA), h = 1/{GA) and c = /{EI) 
where E and G are the Young and shear moduli, A the cross section of the beam 
and / the moment of inertia for fiexion. We used for all simulations the values 
a = 1.0, h = 0.0017 and c = 8.6. When a site is rotated {(pi ^ 0) the beams 
bent accordingly always forming tangentially 90° angles with each other. In this 
way local momenta are taken into account. For a horizontal beam between sites 
i and j one has the longitudinal force acting at site j: Fj = a{xi — xj); the 
shear force: Sj = P{yi — Uj) + ^l{(pi + (fj), and the fiexural torque at site j: 
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Mj = ^l{yi - Uj + l(pj) + 5l'^{(pi - (fj), using a = 1/a, (3 = 1/(6 + c/12) and 
5 = (3{b/c + 1/3). The corresponding equations for vertical beams are similar. In 
mechanical equilibrium the sum over all internal and external forces (torques) act- 
ing on site j must vanish giving rise in the continuum to the Cosserat equations. 
We do not consider here inertial or bulk forces, as for example gravity. 

Before discussing the employed boundary conditions and their physical moti- 
vation it is convenient to describe how we included heterogeneous cohesion prop- 
erties into the fracture model. The concept of a local cohesion strength has been 
used in a number of papers One assumes that a deformed elastic beam con- 
necting sites i and j always breaks above a certain material specific threshold force 
fcoh ('cohesion strength'). If the applied stresses (forces/beam section) are above 
this threshold the beam does break and is removed, i.e. its elastic moduli are 
set to zero. Since the cohesional strength for compression is much higher than 
for tension we assume that compressed beams can never break. If all beams 
have the same cohesion strength the material is homogeneous. Such homogeneous 
states are usually investigated in continuum fracture mechanics using the concept 
of the solid's free surface energy. In the next section we will present some results 
for equal cohesion thresholds. More realistic, however, is the case in which the 
breaking thresholds are distributed randomly according to some probability den- 
sity function, i.e. following a power-law, p{fcoh) ~ floh ^'^^^ fcoh e [0, fmax] and 
r > —1. Negative exponents r arc used to describe strong cohesive disorder while 
large positive exponents correspond to weak disorder. It is convenient to express 
the normalization factor and fmax by the distribution's expectation value {fcoh) 
and the exponent r. We fixed the average cohesion strength for all simulations 
to be {fcoh) = 0.01 and investigated the fracture processes for the exponents: 
r = —0.7, r = —0.5 (strong disorder) and r = -|-oo (no disorder). 

Boundary conditions must be defined on the external edges of the lattice and 
on the internal crack surface. Concerning the external boundary conditions all 

displacements and rotations of the sites on the outer edges are assumed to be 
periodic in horizontal and vertical direction so that the lattice is spanned on a 
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torus. Periodic boundary conditions are preferrable to free boundary conditions 
because we are interested in asymptotic results for large (infinite) systems. We note 
that in this case the system cannot expand globally. The second kind of boundary 
conditions concerns the conditions for forces and torques acting on the inner crack 
surface. While in tensile experiments the crack surface(s) are always stress free, 
in hydrodynamic fracturing these surfaces are loaded by a pressure distribution 
resulting from the invading fluid or gas. In this work we will only consider a 
homogeneous (spatially constant) pressure distribution acting perpendicular along 
the entire inner crack surface. However, in contrast to previous simulations here 
the pressure has strong fluctuations in time. Instead of keeping the inner pressure 
constant during the simulation we consider the case where the fluid flux AV into 
the crack is constant in time. Hence the crack opening volume V which corresponds 
to the total amount of injected incompressibel fluid increases linear in time t, 

V{t) = AVt, AV = const., (1) 

a condition which is close to the situation of industrial hydraulic fracturing. For 
completeness we should mention that the above mentioned equivalence between 
the crack opening and the injected fluid volume does only hold if there is no other 
sink of fluid in the system besides the crack. Although in practice a loss of fluid 
in soils is quite common for seek of simplicity we do not consider here this effect. 

Through eq.(l) the crack volume V increases in time. We will follow the 
evolution of the crack growth under the continous increase of loading in the hole. 
We do this in an iterative way. At each time step the volume increases by a flxed 
amount AV and we estimate from the corresponding elastic solution the stress 
distribution on the crack surface. According to the breaking rule certain beams 
can break at a given time step depending whether their stresses are beyond their 
breaking thresholds or not. When no beam breaks we continue with the next 
time step. However, if beams break they are simultanous and irreversibly removed 
from the set of elastic equations before one proceeds to the next time step. The 
simulation stops when a certain maximum volume Vmax is reached. During the 
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simulation we monitor the acting crack pressure, the number of broken beams and 
the elastic energy of the system. 

In the following we give a more detailed description of the above steps. Our 
simulation starts at time t = 1. At the place into which the incompressible fluid is 
injected (center of the lattice) one vertical beam connecting sites i and j is removed 
from the force and momentum balance equations, i.e. the beam is broken. Since 
we want to simulate the loading of a crack by an injected fluid a pair of opposite 
forces (dipole) of unit strength is applied at the sites i and j pointing into the 
elastic bulk. The (unit) pressure is then just deflned as the acting force per beam 
length I (l = 1). Under the influence of the unit force dipole the lattice becomes 
distorted. Lattice distortions are in general characterized by the displacement 
field u{t) = {xi{t),yi{t)J(^i{t)), which in turn is determined from the boundary 
conditions at time t. In the following we will denote by uo{t) the displacement 
field calculated for a unit pressure Pq. In Ref. we described how to calculate the 
crack volume V from given crack surface displacements u. We denote by Vq (t) the 
crack volume when computed from uo{t). At time t = 1 the injected fluid volume 
V{1) = AV exerts an equilibrium pressure -P(l) = AV/Vo{l) on the crack surface. 
This follows from the linearity of the elastic equations. The elastic solution which 
corresponds to this equilibrium pressure is just u{t) = P{t)uo{t). We note that 
this simple relation only holds if the acting pressure can be considered as spatially 
homogeneous. 

We will consider that only beams along the surface of the inner hole can break. 
In that way only one single crack is generated. At each time step we determine 
for all beams on the crack surface the force f^^ (t) acting along its axis and only if 
this /'"'(t) is positive, i.e. under tension, and larger than the breaking threshold 
fcoh beam is broken, i.e. its elastic constants a, (3, 5 are irreversible set to 
zero. The forces P^{t) must be calculated from the actual displacements u{t). 
Note that at a given time step t either several beams can break simultaneously 
or none at all. The latter is the case if the stresses of all crack surface beams 
are below their thresholds. In such a situation the crack volume for unit pressure 
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does not change, Vo{t + 1) = Vo{t), because the boundary value problem does 
not change, i.e. ■uo(t + 1) = 'Uo(t). Hence it is not nescessary to recalculate the 
equilibrium forces. If a certain number of beams breaks simultanously additional 
unit force dipoles have to be applied at their corresponding neighboring sites i and 
j destroying the previous balance of forces. Then one has to calculate again the 
internal equilibrium of forces. This is done in our case using a conjugate gradient 
algorithm with a precision of e = 10"-'^° (see Eq. (47) in Ref.8). Due to the new 
boundary conditions the elastic solution changes, i.e. Uo{t + 1) 7^ ^0(^)7 ^-nd the 
volume Vci(t + 1) 7^ Vo{t) is recalculated from the new crack surface displacements 
■uo(t + 1), as described in ref. 

At time t+1 the crack pressure takes the value P{t + 1) = V{t + l)/Vo{t+ 1) 
and the elastic solution matching eq.(l) is given by u{t + 1) = P{t + l)uo{t + 1). 
From these displacements one recalculates the forces J*-' (t + 1) for all beams on the 
(new) crack surface, decides which are the next beams to be broken and repeats 
the above described procedure. 

3. Fracture of homogeneous solids 

Although homogeneous cohesive properties in solid systems are rather the 
exception than the rule it is useful to study them mainly because of their relative 
simplicity. Our model should be capable to reproduce some general features of 
simple continuum models for hydraulic fracturing. First we consider the limiting 
case of equal breaking thresholds (fcoh) = fclh ~ ^-^^ ^ lattice of linear size 
L = 150. The crack volume (total amount of injected fluid) is increased at a 
constant rate of AV = 0.05Z^ per time step. The simulations stop when the 
lattices are broken into two pieces. As expected we find the cracks to have a linear 
shape. This is due to the fact that the highest stress enhancements occur at the 
two vertical beams at the tips. Fig. 1 shows the time dependence of the pressure 
P inside the crack in a double logarithmic plot. One can see that on average the 
pressure drops in time and has oscillations on short time scales. At the beginning, 
i.e. for time t = 1, the crack is very small (one vertical broken beam) and one 
needs a high pressure in order to push the fluid of volume AV into the crack. In 
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this particular calculation the two vertical beams on the crack surface are already 
stressed beyond their cohesion threshold. At the next time step, t = 2, these two 
beams are broken and the pressure drops, because the grown crack can now be 
opened much easier than before (see fig. 1). The pressure goes down although 
additional fluid AV^ has been added to the crack at this time step. It is obvious 
that a large crack experiences a much lower pressure than a small crack for the 
same opening volume because the system in that case is locally more stressed. If 
the pressure drops too much (like at time t = 2) the stresses at the two crack 
tips fall below their cohesion value and the crack cannot grow at the next time 
step. By injecting more fluid into the crack the pressure increases linearly in time 
until the cohesion forces can be overcome again. In fig. 1 one sees a sequence 
of this oscillating pressure. In the continuum description a smooth decrease of 
the breaking pressure in time (volume) is expected. Using continuum mechanics 
one can argue that the smallest pressure necessary to extend the crack at a given 
opening volume should behave like Pent ~ V~^/^ in d = 2 Such a relationship 
should only hold for an elastic infinite plate. Because of eq.(l) we can identify 
crack volume with time up to the factor AV. For comparison we have plotted in 
fig. 1 a dotted line having a slope of —1/3. The agreement with our numerical 
values is quite acceptable over one decade. For small and large times we obtain 
pronounced deviations which originate from the lattice structure and the finite 
size of the lattice and from the external boundary conditions. 

It is interesting to consider the temporal evolution of the stored elastic energy 
U. We have calculated the time dependent energy directly by summing up the 
elastic energies of all non broken beams in the system at every time step. In fig. 2 
we show in a log-log plot the time dependence of the elastic energy. Again we 
see an oscillating behavior as discussed above for the pressure. Breaking events 
as for example at time t = 1 and t = 3 decrease the elastic energy while it 
is increased by pushing fluid into the crack. Globally the elastic energy must 
increase because the system becomes more compressed in time. We find that the 
peak energies scale in time as a power law over two and a half decades with an 
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exponent close to 2/3, see fig. 2. This exponent can be understood by calculating 
the work Wcrit — Jq Pcrit{V') dV done by the external forces in order to inject 
the fluid volume V , which yields 2Ucrit = Wcrit oc V'^^'^ . This agrees well with our 
numerical findings. Again finite size effects lead to deviations from the behavior 
of an infinite continuum. 

4. Cracks in heterogeneous solids 

In the following we will consider the influence of strongly heterogeneous cohe- 
sive strengths on the hydraulic fracturing process. It is well known that fracture 
processes in disordered solid systems show a rich phenomenology concerning the 
crack geometry as for example crack branching and crack deflectiont^l . There exist 
phenomena, as for example the appearence of microcracks in solids, which cannot 
be explained by equilibrium thermodynamics without considering the heterogenity 
of physical properties. Experimental and theoretical investigations during the last 
decades have manifested that the overwhelming number of fracture phenomena in 
solid systems is strongly influenced or even controlled by inhomogeinities. 

We will investigate the statistical properties of breaking sequences during 
hydraulic fracturing and their correlations in space, time and magnitude. Similar 
quantities are frequently considered in the analysis of earthquake occurrences 1^°'^^] 
and of acoustic emission (AE) during microfracturing of rocks '■'^^1 or of technical 
materials t^^^ . Acoustic emission records for the hydraulic fracturing in geothermal 
wells have been published for example in Ref. ^^^1. 

In our simulations we have considered threshold distributions for two different 
exponents r = —0.7 and r = —0.5 (see the previous section for the definition of the 
distribution). Fig. 3 shows a typical hydraulic crack pattern grown in a medium 
with strong disorder (r = —0.7) on a lattice of linear size L = 150. The crack con- 
sists of 629 broken beams after 1500 time steps. Apparently smaller crack branches 
appear along larger branches, a geometrical property which is typically observed 
for self similar (fractal) structures. We have evaluated the fractal dimension(s) of 
the generated hydraulic cracks by considering the relationship between the typical 
crack radius R{N) and the number of broken beams N. For this we calculate the 
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squared radius of gyration R'^{N) = ~ ^o)'^ with ro = J2i '^i- Finally 

the radius R{N) is averaged over a number of independent configurations in order 
to obtain results for typical cracks. In fig. 4 we show in a log-log plot the number 
of broken beams N versus the typical crack radius R{N) for the two statistics of 
the heterogenities characterized by exponents r = —0.7 and r = —0.5. We find 
in both cases power laws N oc R^f with df > 1, giving evidence for fractal crack 
growth. The exponent dj is called the fractal dimension of the crack and takes for 
r = -0.7 (o) and r = -0.5 (+) the values df = 1.44 ± 0.10 and df = 1.39 ± 0.10 
respectively. These values for df are consistent with the fractal dimension found 
in the two dimensional experiments on hydraulic fracturing of viscoelastic clays . 
One detects, however, from both curves in fig. 4 a crossover at large i? to a lower 
slope df ^ 1.25. It is interesting to note that the crossover radii i?x depend on 
the width of the threshold distribution. The cracks for stronger disorder r = —0.7 
show a higher crossover radius i?x ~ 18 than those for r = —0.5 with R^ ~ 10. 
This is qualitatively in agreement with the observation that broader threshold 
distributions have larger homogenization volumes I"*^^]. In our case it is, however, 
likely that the crossover behaviour is an effect of the finite size of the lattice. It 
has been argued ^^^^^^^ that for logarithmically diverging threshold distributions, 
i.e. for r = — 1, a fracture process essentially reduces to a percolation problem. 
We note that this is not necessarily the case here because the employed breaking 
law is asymmetric with respect to tension and compression. 

5. Bursts and temporal correlations 

Since the fractal nature is a finger print for infinite-range correlations in the 
crack geometry, one can also ask how the breaking events are correlated in time. 
To illustrate this question we show in fig. 5 the complete breaking sequence of the 
crack displayed in fig. 3. We have plotted the number of beams broken between 
two consecutive time steps as a function of time. Most striking is the fact that 
the breaking process is very discontinuous in time. There are large time intervals 
in which no breaking occurs at all. During such time intervals of quiescence all 
beams on the crack surface are stressed below their cohesion thresholds and the 
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acting pressure increases linearly in time. The second striking fact is that if a 
breaking event happens after a period of quiescence it usually triggers a sequence 
of consecutive breaking events (temporal clustering). We will call such sequences 
bursts. The bursts themselves are, as one can see in fig. 5, unequally distributed 
in time. They occur relatively often for small times and become more rare later. 
However, while the bursts are narrow for early time steps they typically become 
broader with increasing time. The numbers of simultanously broken beams per 
time step also exhibit particularities. Let us call these numbers the magnitudes of 
the breaking events. Broad bursts typically consist of few events of high magnitude 
and much more of low magnitude. Fig. 6 shows the magnifications of two bursts 
from fig. 5. The ordinate gives the number of simultanously broken beams during 
a time step and the corresponding 'released energies' as well. The definition of 
the released energies is given further below in Sec. 6. We see that the peak energy 
releases do not temporarly coincidence with the peak numbers of broken beams. 
In fig. 7 we show the time dependent pressure belonging to the breaking process 
shown in fig. 5 and fig. 3 respectively. 

The reader familar with magnitude records of earthquakes or of accoustic 
emission records from laboratory experiments will recognize some ressemblance 
with our data. 

In the following we will investigate more closely the temporal clustering of 
breaking events. We have calculated the lifetime r for each burst as the time that 
elapses between the first and the last breaking event of a burst. Fig. 8 shows the 
(unnormalized) histogram of burst lifetime in a double logarithmic plot. Small 
bursts occur relatively often while larger bursts are less frequent. The statistics 
are made over 729 bursts from 60 samples for r = —0.7 and over 862 bursts from 
53 crack simulations for r = —0.5. The largest detected burst has a lifetime of 
approximately r = 140. We note that for large lifetimes, r > 30, the statistic 
becomes unreliable because the occurrences, n(T), become too small, n(T) < 10. 
This is also due to the fact that all simulations were stopped after tmax = 1500 
time steps. In order to extract more information from the lifetime distribution we 
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consider in fig. 9 the less noisy cummulative probabifity distribution p{t) that a 
given burst has a hfetime shorter or equal to r. With an intermediate range we 
observe in fig. 9 that the cummulative probability p(r) = X]i=i ^(0/ X]i=i^ 
seems to obey a power law, p(r) oc t^~'^, with rj — 0.54 ± 0.15. Hence the lifetime 
distribution of bursts also follows a power law in this regime, 

n(T) cx T"''. (2) 

It is remarkable that the two curves in fig. 9 give nearly the same exponent rj 
although the exponents r from their threshold distribution are diff'erent. This 
indicates that the underlaying generating mechanism for the temporal clustering 
of breaking events might be universal. In fig. 9 one also sees strong deviations 
from the power law for small (r < 4) and for large (r > 30) bursts particularly in 
the case of r = —0.7. The existence of an upper cut off is plausible, because we 
stopped each crack growth after tmax = 1500 time steps. This artificially lowers 
the number of large bursts. The lower cut off is quite common in cluster statistics, 
known as 'corrections to scaling'. 

Next we consider the lifetime distribution q{T) of quiet intervals. This dis- 
tribution characterizes the arrangement of bursts on the natural time scale t and 
has attracted interest in seismology because of its role for possible predictions of 
earthquakes t^^l. It has been observed that if certain regions undergo high magni- 
tude earthquakes they often show thereafter rather long periods of seismic quies- 
cence (inactivity). We have investigated this statistics from our hydraulic fractur- 
ing calculations. In fig. 10 we consider in a semi-log plot the cumulative probability 
Q{t) that two consecutive bursts are separated by a time interval of quiescence 
shorter or equal equal than r. More precisely, we define the width of the interval 
T by the number of breaking free time steps between two consecutive bursts. For 
intermediate time scales the probability follows, as shown in fig. 10, a logarithmic 
dependence, Q{t) oc Inr. Thus the probability q{T) to find a quiet time interval 
of width T between two adjacent bursts scales in this regime as, 

g(r) « ^. (3) 
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Interestingly both distributions n(T) and q{T) scale for intermediate times as power 
laws, however, with different exponents. 

So far we have considered the lifetime statistics of bursts and of quiet inter- 
vals. Valuable information about the temporal evolution of hydraulic cracks can 
be gained by considering the time correlations between the breaking events. A 
priori it is not obvious which of the breaking events are causually connected. In 
our model the elementary consecutive breaking events define the shortest accessi- 
ble time scale. From acoustic emission experiments and seismic records, however, 
the existence of a background noise is well established. Events are usually only 
considered if their magnitude exceeds the background noise by orders of magni- 
tude [^^'^^]. The cut off of the background noise leads to discrete magnitude-time 
sequences. This is similar in our model because the discrete nature of the beams 
put a lower cut off on the possible size of an event. Experimentally, the variations 
in emission magnitudes are so strong and the line- widths are so small that one can 
consider acoustic emission 'events' as delta peaks. We note that in experiments 
the significant breaking events are already temporal clusters (bursts). Hence one 
can investigate the time correlations on the scale of a typical line-width (burst 
lifetimes) or on a larger time scale. First we will investigate the case in which only 
correlations within the bursts are considered. We have calculated from all bursts 
the probability distribution 6(t) of time intervals t = \ti — tj\ between all possi- 
ble pairs of breaking events belonging to the same burst. We show in fig. 11 in a 
semi- logarithmic plot the distribution of time intervals r. One clearly sees an expo- 
nentially decreasing probability (two-point correlation) in time, 6(t) oc exp (— ar). 
In this plot one does not find any power-law behavior for small times, r < 30, as 
one might expect from the power-law scaling of burst lifetimes, see eq.(2). In our 
case this is due to the fact that the 'large' bursts (r > 30) having a bad statis- 
tics completely dominate the correlations also for small r. This appears to be an 
additional complication resulting from the very broad lifetime distribution n(T) of 
bursts, i.e. the small value for the exponent i] ~ 0.5. We can, however, as dis- 
cussed above consider the correlations between all breaking events regardless of the 
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bursts they belong to. This is conveniently done by calculating the histogram 
of time intervals r between all possible pairs of breaking events r = \ti — tj\ from 
a given time record. In fig. 12 we show in a double logarithmic representation 
the normalized histogram of these time intervals. We find that the probability 
g{T) to detect two breaking events seperated by a time interval r decreases on 
intermediate scales (15 < r < 130 for r = —0.7) as a power law, 

g{T)(XT-'', At = 0.94 ± 0.20. (4) 

This corresponds to Omori's well known law of aftershocks which has been 
first formulated in 1894 for earthquakes It says that the probability for an 
aftershock following a large earthquake decays as l/r'* in time with an exponent 
K close to one. Omori's law has been verified from earthquake catalogs for after- 
shocks series ranging from a few hours to a couple of years after the main event. 
The empirical value for the exponent k, lies between 1.0 and 1.4 1^°]. Though the ap- 
pearance of earthquakes and aftershock series underlies quite different constraints 
concerning the rheology and the crack opening mode than in our two dimensional 
model we find a value for k very close to the empirical one in three dimensions. 
Recently it has been argued [^^^ that the value n = 1 should hold quite generally 
for nonlinear fracture problems, such as earthquakes, independent on the spatial 
dimension 1^°] . Hence one might expect the basic mechanism for bursts sequences, 
respectively aftershocks, to be universal due to self organisation of rupture. Phys- 
ically one might argue that a main breaking event with high magnitude triggers 
smaller magnitude events which in turn create even smaller events ad infinitum if 
the cohesion properties are sufficient heterogeneous. However, in our case during 
a single burst the solid undergoes a quite complicated stress redistribution pro- 
cess until all microstructural elements (beams) are below their cohesion threshold. 
When the crack grows into new regions it may be stopped due to a region with 
particular high breaking thresholds. Then the burst terminates. As the crack 
becomes longer its critical opening volume necessary to overcome the thresholds 
increases even more because larger cracks can be opened easier than small cracks. 
This explains why the intervals of quiescence become longer for larger cracks. 
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6. Released energies, acoustic emission 

So far we have considered the spatial and temporal correlations of breaking 
events during hydraulic fracturing. We have found qualitative and quantitative 
similarities to the seismic clustering of earthquakes and to acoustic emission labo- 
ratory experiments. However, we have so far not discussed the correlations of the 
magnitudes (intensities) of breaking events. For this the magnitude occurrence 
relationship is of central importance t^^^. It has been found empirically that the 
cumulative occurrence H{m) of earthquakes of magnitude larger than m follows 
the celebrated Gutenberg- Richter law, H{m) oc 10"*"" l^^l The 'b-value' found 
in Gutenberg- Richter's law is often close to one t^*^!. The magnitudes of seismic 
events are defined by the logarithmic wave amplitudes. The energy release, 5U , 
is usually considered to be proportional to the squared amplitude of the earth- 
quake. One typically obtains a cumulative occurrence-energy relation of the form, 
H{dU) oc SU~^'^, expressing the presence of earthquakes on all energetical scales 
[11,40] Yias been shown recently [^"^^ that acoustic emission (AE) in the ultrasonic 
frequency range due to microfracturing of synthetic plaster also exhibits a power 
law for the cumulative occurrences with an exponent close to the one found in 
the Gutenberg-Richter law. Similar findings from AE records have been published 
earlier for the microfracturing of rocks I^^^ . 

In the following we will discuss some energetical aspects of our model. The 
thermodynamical description of crack spreading phenomena in heterogeneous en- 
vironments becomes quite complicated and we will calculate the involved energies 
exclusively from the mechanical equilibrium conditions and not from thermody- 
namical considerations. In order to determine in our model the amount of energy 
released by a given breaking event one has to calculate the stored elastic energy 
U (t) at each time step. This energy is, as already mentioned, determined prior 
to breaking from the sum over the energies of all non broken beams. We define 
the rate of change of elastic energy as 5U = —dU/dt = —{U{t + 1) — U{t)). If 
the rate 5U is negative then the lattice has increased its elastic energy. We have 
verified that the overwhelming number of breaking events are detectable by the 
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condition SU > 0. There are very few events that show a negative energy flux 
because the value of 5U describes the sum of two effects. First, it describes the 
lowering of elastic energy due to a breaking event. Second, it describes the increase 
of elastic energy due to the added amount AV of fluid per time step. However, 
if the cracks and their corresponding opening volumes become large compared to 
AV the fraction of non-detected events becomes very small. Experimentally the 
energy dissipation due to crack growth has at least two contributions. On one 
hand the creation of additional crack surface consumes energy. On the other hand 
elastic waves are emitted from the crack tip(s) (acoustic emission (AE)). Because 
we do not consider the full dynamic equations it is a priori not clear which fraction 
of the dissipated energy to assign to the 'acoustic emission'. We propose that, as a 
flrst approximation, the AE energy per event is proportional to the rate of change 
of elastic energy dU of the system. We show in fig. 13 in a semi-log plot the cumu- 
lative occurrence H{5U) for breaking events of energetical magnitude larger than 
5U . We find an exponential relationship, 

oc exp(-/35C/). (5) 

The flat tail for large 5U is due to low statistics. This result is at variance with what 
is expected from the Gutenberg-Richter law and AE measurements in laboratory 
experiments. A possible source for this discrepancy could be our deflnition of the 
released energy: We proposed that the AE be proportional to the total energy 
relaxation. It could, however, be that the fraction of energy emitted acoustically 
is a more complicated function of the energy. Other sources for the deviations 
between our model and experiments could be the periodic boundary condition 
and the two dimensionality of the lattice used in our simulations. 

7. Conclusions 

We have presented a lattice model for hydraulic fracture which takes into ac- 
count the particular boundary conditions inside the hole and the quenched nature 
of the heterogeneities in the rock. We drive the process by increasing the crack vol- 
ume linearly in time (constant flux). The pressure fluctuates erratically similarly 
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to what is measured at boreholes in the field. The crack shapes are fractal and the 
fractal dimension agrees well with measurements performed in Orleans with Hele 
Shaw cells. The sequence of breaking events is organized in bursts which have a 
life time distribution and a distribution of quiescence times that are power laws. 
This indicates that self-organized criticality (SOC) ^^^^ takes place. Events in- 
side a burst seem uncorrelated while long range correlations between bursts follow 
Omori's law for aftershocks. The distribution of released energies does not follow 
a power-law within the numerical accuracy of our calculations, shading doubts on 
the simple hypothesis that acoustic emission signals are simply proportional to the 
released potential energy. 

This work is still rather preliminary if the full reality of hydraulic fracturing 
in oil or geothermal reservoirs is to be described. Real soils are three dimensional 
and the distribution of heterogeneities usually follows a WeibuU distribution with 
a material dependent exponent. The pressure of the fluid in the crack is not 
constant but depends on the distance from the injection site and the geometry 
of the crack. The system is nearly infinite in size and the restrictions in the 
total volume of the system as imposed by the periodic boundary conditions are 
not realistic. Still many features have been found in this paper that agree with 
experimental measurements and we think that including more details into the 
model, which is rather straightforward, could help simulate hydraulic fracturing 
rather accurately in the future. 
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Figure 1: Log-log plot of the pressure P inside the crack versus time t for 
homogeneous cohesion {r — > +oo) of strength (fcoh) = 0.01. Points at subse- 
quent time steps are connected by straight lines. The dotted line corresponds 
to a slope of —1/3 as predicted by continuum mechanics. The linear lattice 
size is L = 150. 
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Figure 2: Log-log plot of the stored lattice energy U versus time t for homo- 
geneous cohesion {r — > +oo) of strength (fcoh) = 0.01. Points at subsequent 
time steps are connected by straight lines. The dotted line is a guide to the 
eye of slope 2/3 as predicted by the continuum theory. The lattice size is 
L = 150. 



Figure 3: Typical hydraulic crack for strong disorder, r = —0.7 and {fcoh) = 
0.01, on a lattice of size L = 150. The crack consists of 629 broken beams 
after 1500 time steps using a volume increment AV = 0.05. The injection 
point is the center of the lattice. We have used periodic boundary conditions 
in vertical and horizontal directions. 
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Figure 4: Log-log plot of the number N of broken beams versus the radius 
of gyration -R(A^) for two different distributions of breaking thresholds: (o) 
r = —0.7, averages over 60 cracks, fractal dimension df = 1.44 ± 0.10; (+) 
r = —0.5, averages over 53 cracks, fractal dimension df = 1.39 ±0.10. For all 
simulations we have used a mean cohesion value (fcoh) = 0.01 and a linear 
lattice size L = 150. 
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Figure 5: Record of the breaking sequence in time corresponding to the crack 
displayed in fig. 3. In this plot we have defined the magnitude m(t) as the 
number of simultanously broken beams at a given unit time interval. The 
simulation stopped after 1500 time steps. Note the temporal clustering of 
breaking events and the large time intervals of quiescence. 
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Figure 6: Magnification of two bursts from fig. 5, (a) at time t ^ 650 and 
(b) at time t ~ 1210. We have plotted the number m{t) of simultanously 
broken beams (thick line) as well as the corresponding released energies SU 
(thin line) (compare Sec. 6 for the definition oiSU). The energies have been 
scaled by a constant factor. 
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Figure 7: Linear plot of the pressure P inside the crack versus time t. The 
pressure was obtained from the simulation corresponding to fig. 5. 
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Figure 8: Histogram for the burst lifetime r in a double logarithmic repre- 
sentation. The occurrences nij) for bursts of size r have been calculated for 
r = —0.7 (o) and for r = —0.5 (+). The largest detected burst is of size 
T = 140. The simulations have been stopped after tmax = 1500 time steps. 
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Figure 9: Log-log plot of the cumulative probability distribution pij) that a 
burst selected at random is shorter than r. The distribution follows a power 
law, i.e. pij) oc t^~^ with rj = 0.54 ± 0.15 for intermediate burst lifetimes. 
Symbols and parameters as in fig. 8. 
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Figure 10: Semi-log plot for the cumulative probability Q(t) to find between 
two adjacent bursts an emission (breaking) free time interval r. Note the 
logarithmic dependence, Q{t) oc Inr, i.e. the probability qij) to find a quiet 
interval of size r scales as qij) = d^Q oc 1/r. Symbols and parameters as in 
fig- 8. 
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Figure 11: Probability distribution 6(t) to find two breaking events being r 
time steps apart and belonging to the same burst. From the semi-log plot one 
sees exponentially decreasing correlations, 6(t) cx exp(— ar). Symbols and 
parameters as in fig. 8. 
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Figure 12: Probability distribution (/(r) to find two breaking events to be 
T time steps apart. The distribution follows for intermediate time scales a 
power law (Omori's law), g{r) cx r"*^ with the value k = 0.94±0.20. Symbols 
and parameters are the same as in fig. 8. 
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Figure 13: Semi-logarithmic plot for the cumulative occurrence H{6U) that 
a given breaking event is of energetical magnitude larger than 5U. We find 
for the cumulative occurrence an exponential decay, H{5U) cx exp{—f35U). 
Symbols and parameters as in fig. 8. 



